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Abstract. The Baryon Acoustic Oscillations (BAO) in the large-scale structure of the uni- 
verse leave a distinct peak in the two-point correlation function of the matter distribution. 
That acoustic peak is smeared and shifted by bulk flows and non-linear evolution. However, 
it has been shown that it is still possible to sharpen the peak and remove its shift by undoing 
the effects of the bulk flows. We propose an improvement to the standard acoustic peak 
reconstruction. Contrary to the standard approach, the new scheme has no free parame- 
ters, treats the large-scale modes consistently, and uses optimal filters to extract the BAO 
information. At redshift of zero, the reconstructed linear matter power spectrum leads to 
a markedly improved sharpening of the reconstructed acoustic peak compared to standard 
reconstruction. 



1 Introduction 

The Baryon Acoustic Oscillations (BAO) are an important target for studies of large-scale 
structure (LSS) in the universe. They appear as an acoustic peak in the 2-point correlation 
function of the galaxy distribution [1]. The scale of the acoustic peak corresponds to the 
distance sound waves have traveled in the photon-baryon fluid until recombination. That 
scale (~ 150 Mpc) is well in the linear regime, and bulk matter flows and non-linear physics 
can affect the scale only at a percent level [2]. Yet, these effects introduce a ~ 10 Mpc 
broadening of the acoustic peak, which in turn degrades the accuracy in the determination 
of the acoustic scale. Most of this smearing is due to large-scale flows [3], whose effects can 
be reversed at least partially [4-6] in a process known as reconstruction. As LSS surveys 
become larger and more complete, increasingly more advanced algorithms to compensate 
for the degradation of the acoustic peak will be needed to fully utilize the observations. 
Ongoing and future surveys of LSS are expected to achieve a sub-percent level accuracy in 
the determination of the acoustic scale, making the development of such algorithms the more 
urgent. 

In this paper we propose a new scheme for the reconstruction of the acoustic peak, 
the goal of which is to recover the linear density field starting from the non-linear matter 
density field. The method relies on a recent model [7] for the non-linear density field, where 
one approximates the trajectories of Cold Dark Matter (CDM) particles as given by a scale- 
and time-dependent rescaling of the particle displacements in the Zel'dovich approximation 
(ZA) [8] . That model allows for the construction of optimal filters, with which to iteratively 
process the matter density field, so as to compensate for the smearing and shifting of the 
acoustic peak induced by the large-scale matter flows. This is done without introducing any 
free parameters, unlike the standard reconstruction scheme which relies on finding an optimal 
smoothing scale by trial-and-error. 

An important ingredient of the method described in this paper is its treating of the 
large-scale modes. Standard reconstruction recovers them in a Zel'dovich-like approximation, 
which implies (as we will see below) that one can choose to recover well either the large-scale 
modes (k < 0.2/i/Mpc for redshift of z = 0); or the short-scale modes (k > 0.2/i/Mpc for 
z = 0) by degrading the large-scale reconstruction. Our method ameliorates this trade-off 
considerably. This results in a significantly improved sharpening of the acoustic peak in the 
matter correlation function. 

In this paper we apply our reconstruction scheme only on the matter distribution in 
real space, as we rely on previous results [7] (see also [3]) to model the CDM matter density 
field. Once those results are extended to cover biased tracers, it should be straightforward to 
write down the proposed reconstruction scheme for biased tracers of the underlying density 
field. However, it is not yet clear how much the reconstruction will be improved in that case. 

In Section 2 we give a review of standard reconstruction, while in Section 3 we give a very 
brief overview of our method. Then in Section 4 we derive a pseudo-optimal 1 reconstruction 
algorithm, given the model of [7]. We find that that algorithm converges slowly, and therefore 
in Section 5 we present an efficient algorithm which provides a good initial guess for the 
pseudo-optimal algorithm in Section 4. We compare the proposed reconstruction scheme 



1 By "pseudo-optimal" we mean that our derivation would have produced an optimal estimator if the 
fields involved were Gaussian. However, the density field in the (mildly) non-linear regime is non-Gaussian. 
Therefore, all we can say for sure is that our method provides better reconstruction results than the standard 
reconstruction scheme - a result which we obtained only after actually applying the method to N-body 
simulations. 



with standard reconstruction in Section 6, and then we summarize our results in Section 7. 
In the Appendices we give the details behind our numerical implementation of the proposed 
reconstruction scheme. 

2 Standard reconstruction 

Let us start with a review of standard reconstruction [4-6], which will allow us to set up 
some notation. 2 The exact position, x p , of particle p in Eulerian space (obtained e.g. from 
N-body simulations) is given by 

x P (q,v) =Q P + s(q P ,ri) , (2.1) 

where q p is the Lagrangian particle position; r] is conformal time; and s is a stochastic vector 
field. As an example, at linear order in the overdensity we have 

8 q ■ s z (q, n) = -S L (q, r,) = -6 L (q, r, )D(r,) , (2.2) 

where 5^ is the linear density field, D is the growth factor such that D(r/o) = 1. The above 
equation corresponds to the Zeldovich approximation (ZA) [8], hence the subscript z. 

The overdensity field, 5[aj p ](a;), corresponding to the distribution of particles x p , is a 
functional of the particle positions and is proportional to the sum over the number of particles 
falling in a volume element around x. As an illustration, starting from a set of CDM particles 
sitting at coordinates x p given in Eulerian space (a;), one could easily obtain the overdensity 
(5) using for example a Cloud-in-Cell (CiC) assignment. In the continuum limit, this particle 
representation results in 

8[q p + s(q p , n)] (x) = / d 3 qd D (x - q - s(q, r?)) - 1 , (2.3) 

where 5d is the Dirac delta function. 

Standard reconstruction then works as follows. The density field, 5, is Gaussian smoothed 
(with a filter Wq = exp(— £ 2 fc 2 /2)), and then a large-scale displacement field, sr(x), is ob- 
tained to linear order from eq. (2.2): 

s R (k)=i^W G (k)5(k), (2.4) 

where R stands for "Reconstructed", and the wavevectors above correspond to Eulerian 
Fourier space. One then displaces the particle positions, (2.1), by — sr, obtaining the follow- 
ing density field: 

S s = 6[x p - s R (Xp)} , (2.5) 

which to leading order reduces to (1 — Wg)5l- Note that x p — sr{x p ) is an approximation to 
the Lagrangian positions of the particles, since we have subtracted the (dominant) large-scale 
displacements. Therefore, 5 S is really an approximation to the short-scale part of the linear 
overdensity in Lagrangian space, 5i{q). 



2 In this paper we compare only against the one-step ZA procedure described in [5]. Those authors investi- 
gate the effect of incorporating second-order Lagrangian perturbation theory in their reconstruction scheme, 
as well as performing several iterations of their method. The improvements they find are only marginal, and 
therefore we do not focus on those modifications here. 



To get the large-scale modes in the standard reconstruction method, one moves a uni- 
form random particle field, q u (u running over the random particles), by sr to obtain 

<^A = S[q u + s R (q u )] , (2.6) 

which to leading order reduces to Wq$l- The resulting approximation, 6^, to the linear 
density is then given by 

Sf = 5 S + 6 A , (2.7) 

where "st" stands for standard reconstruction. 3 

Note that eq. (2.6) effectively mixes Lagrangian and Eulerian coordinates. This is 
because q u can be treated as particle positions in Lagrangian coordinates, since S[q u ] is 
uniform by construction. Yet, sr is an Eulerian vector field, which we evaluate at the 
Lagrangian instead of the Eulerian positions. This difference of course is irrelevant at leading 
order, but we are interested in the higher-order corrections. Indeed, if we cared only about the 
linear order, we could have simply written 5l = 5, which can hardly be called a reconstruction. 

Even if sr was the correct Lagrangian space displacement field, what eq. (2.6) would give 
us is the density field in the ZA corresponding to the displacements sr. Therefore, the cross- 
correlation coefficient between 5l and 8\ exhibits a Gaussian decay on a scale approximately 
given by the large-scale rms particle displacements. One finds the same behavior for the 
cross-correlation coefficient between 5l and 5 (see e.g. [3]). It is the presence of Wq that 
makes the former decay smaller by reducing the rms particle displacements, which is why 
standard reconstruction works at all. 

3 Brief overview of our method 

From our discussion of standard reconstruction so far, we can see that there are two main 
places where standard reconstruction can be improved: 1) One can try to obtain optimal 
filters for processing the non-linear density field, instead of relying on filters chosen by hand 
(Wg)\ 2) One can treat the large-scale modes consistently by not mixing Eulerian and La- 
grangian coordinates, and by obtaining the large-scale density field as minus the divergence of 
the Zel'dovich displacements, instead of by using the Zel'dovich-like approximation of (2.6). 
These issues are solved in our reconstruction scheme outlined below. 

Before we proceed with the description of our method, one should note that there are two 
main approaches to reconstruction: "forward" and "backward" . In a forward reconstruction, 
we try to match the non- linear 5 with a 5[s z ] corresponding to our reconstructed s z - a fitting 
which we perform in Eulerian space. In a backward reconstruction, we try to move around 
the data particles, x p , responsible for the non-linear 5 in such a way so as to get as uniform 
a field, S[x p — s z ], as possible. In other words, in a backward reconstruction we try to match 
the uniform initial conditions, and therefore this kind of fitting can be thought of as being 
performed in Lagrangian space. 

One should note that if one performs backward reconstruction on real galaxy catalogs, 
one implicitly assumes that the data particles (the galaxies) cover Lagrangian space uniformly 
- a highly non-trivial assumption for biased tracers. So, a forward scheme should be the 

3 In the actual standard reconstruction, the plus in both (2.6) and (2.7) is replaced with a minus. This is 
because in the standard reconstruction scheme one displaces the random points by —sr to get 8a- 



scheme of choice for real-life data. Yet, in experimenting with different algorithms, we found 
that backward schemes are generally both more stable and faster converging. 

Therefore, the reconstruction scheme proposed in this paper consists of two stages. 
In the first stage we use a well-motivated but yet ad hoc backward scheme for efficiently 
constructing a good guess for s z . That guess is then used as an initialization for our forward 
quasi-optimal second stage of the reconstruction. 

In the first stage (Section 5) we solve iteratively for a Lagrangian displacement field 
d(q) which moves the data particles (e.g. galaxies) to form a uniform distribution (see 
(5.1)). We then relate the displacement field d to the ZePdovich displacement, s z , using 
a linear model (see (5.4)). Therefore, the large-scale s z is given by a Wiener-filtered d, 
while the large-scale reconstructed linear 5 is given by minus the divergence of the large-scale 
s z . This is in sharp contrast with standard reconstruction, where the large-scale modes are 
reconstructed in a ZePdovich-like approximation (we will return to this point later on). The 
Wiener filtering, however, loses information about the short scales, which we then recover in 
a manner analogous to standard reconstruction, (2.5). 

The resulting reconstructed s z (given in (5.9)) is then used as an initial guess for the 
quasi-optimal forward second stage (Section 4) of our reconstruction scheme. In that stage, 
we try to maximize the likelihood function given by (4.3), which uses information about the 
statistics of the initial density field, and tries to match the non-linear 5 to an LPT-inspired 
model for the non-linear density field. The result of that procedure is an s' z (see (4.12)), 
which is very well-correlated to the true Zel'dovich displacement field. Yet, s' z is a non-linear 
Wiener-like filtered displacement field, and thus it is suppressed at short scales. We recover 
them by filtering s' z with an inverse Gaussian (see (4.14)) to obtain s z . Our final result for 
the reconstructed linear density field is then given by minus the divergence of s z , according 
to (2.2). 

In the next two sections we will describe in detail the two stages of our reconstruction 
scheme. First we start our discussion by deriving the quasi-optimal second stage of our 
method. The second stage is independent on the first stage of the reconstruction, which only 
supplies an initialization for the second stage. Then, in Section 5 we give the detailed recipe 
for our choice of initialization, constituting the first stage of the method. 

4 Improving standard reconstruction 

In this section, we will derive a pseudo-optimal iterative scheme for performing the BAO 
reconstruction. We will start by writing down a model for the CDM density field, presented 
first in [7] . Then, we will write down the corresponding Likelihood function, £, and from there 
we will obtain an iterative algorithm, with which we will solve for s z , which maximizes C 
That iterative algorithm will be our reconstruction scheme, because knowing s z is equivalent 
to knowing 5l according to (2.2). 

4.1 Modeling the density field in the mildly non-linear regime 

To make progress we use the following split for the non- linear density field proposed in [7]: 

5{k) = Rs(k)Sx(k) + S MC (k) , with Si = 5[x p = q p + (R z * s z )(q p )} , (4.1) 

where 5mc is a mode-coupling term, such that (5i(k)5^ 4C (k)) = 0. Here k is in Eulerian 
Fourier space and the * denotes convolution in Lagrangian space. We introduced Rg(k) and 



i? z (k) - the so-called density and displacement transfer functions, which are given in [7]. 
Here k denotes a wavevector corresponding to Eulerian Fourier space, while k denotes a 
wavevector in Lagrangian Fourier space. These transfer functions are given by a linear fit to 
the density and the displacement fields respectively: 

Rsik) - mm ' zi) ~ d^(k)i 2 ) ■ ( } 

These transfer functions can be cheaply obtained by calibration with simulations because of 
their small sample variance. We refer the reader to [3] and [7] for further discussion on this 
point, as well as to [7] for the details of the above model for 5 used in this paper. 

4.2 The Likelihood Function 

The likelihood function, C, for s z given the model above is: 

where we made the difference between Lagrangian and Eulerian wavevectors explicit. The 
power spectrum Pmc is the power spectrum of 5mc-> but can include other noise terms, such 
as Poisson noise if 8 is sampled with discrete tracers. In writing C above, we made the 
simplifying assumption that 5mc is Gaussian. 

The pseudo-optimal 4 s z for the Likelihood function above can be found by taking the 
functional derivative of log C with respect to s z and setting to zero. For that we need the 
following derivative 



d5{{k) f d 3 



f dVjk\ q+ {R**s*)( q J) ik . e -*-«7^( k ) 

J (2-kY 



(4.4) 



where we denote functional derivatives by d so that we can distinguish them from the over- 
density Here again the Lagrangian wavevector is denoted with k, while the Eulerian - by 
k. To derive the above equation we wrote down S±(k) by Fourier transforming (2.3) with 
s — > Rz * s z . Thus, for a quantity A(k) we can write: 

f d ^ k ^M A{k) = f ^^A(k)e ik i^ R ^^)ike-^R z (k) 

J os* z {k.) J (2tt) 

= J 7^3 d 3 x5 D {x -q-(R z * s z )(q))V x A(x)R z (k)e-^ . (4.5) 

Combining (4.5) and (4.3) we obtain that the maximum of the likelihood function is at s z 
that satisfies: 

l dlog£ , . 

= "2^(k) = (46) 



,2 



s z (k)-R z (k) I -fi- e - ik «V 



P L (k) " ' ^ ' J (27T) 



d A ke 



PMc(k) 



4 We have assumed a Gaussian 5uc, which makes our calculation no longer optimal. 



where x is evaluated at x = q + (R z * s z )(q), and one should note that 5\ is a functional of 
s z . The integral in q above gives the Lagrangian Fourier transform of the Eulerian gradient 
of the quantity in the square brackets. Thus, we can write the above expression as 



l s z (k) = R z (k)P L (k)T kq T-£ 



ik 



R 5 (k)(6(k) - R 5 (k)6 1 (k)) 
PMc(k) 



(4.7) 



where J- kq denotes the Fourier transform from q to k, while J-~ k denotes the inverse Fourier 
transform from k to x. The reason J"kq^xk does not equal 1, is because between those two 
transforms we map the coordinates according to x = q + (R z *s z )(q). At first order, however, 
the difference between x and q (and between k and k) is irrelevant and 5\ = — V • {R z * s z ), 
and thus we obtain: 



ik • s z (k) = 



(R z Rs) 2 P L 6(k) 

Pmc + {RzR&) 2 Pl X R z Rs 



(4.8) 



which is the standard Wiener filter result for a linear model relating the displacement and 
the non-linear overdensity. Thus, we should think about the result for s z in (4.7) as coming 
from a non-linear Wiener-like filtering of 5. Note, however, that in what follows we use the 
exact result given by (4.7), and not the approximation given by (4.8). 

4.3 Maximizing the Likelihood function iteratively 

We would like to find the maximum of the Likelihood function iteratively using the usual 
Newton-Raphson method. Thus, we write the updated value s' as: 



4i(k) 



m 



d 2 log C 



6 D (o) \ds z (k)ds* z (k) 



-i 



dlogjC 

■4 9K~(F) 



(4.9) 



where we used homogeneity (which is why we have the factor of 1/<5d(0) = (2vr) 3 /l/, where 
is the zero wavevector, and V is the (Lagrangian) volume of the survey/N-body simulation); 
and as usual we used the Fisher matrix instead of the Hessian of the Likelihood function in 
order to simplify the inversion. 

Next we need take the second derivative of the likelihood function in analogy with the 
first derivative. We find: 



1 



d 2 log £ 



28 D (0) ds g>i (k)dtf Z)j (k) 

R 2 (k 



Pl( 



+ 



(4.10) 



, d 3 A:d 3 ^ 3 ^27r)- 6 ^ ( ^ T A : ^ J e ife -(' ?+ ^*^)( ,? )-^-^*^)^) e - lk -(' ? -^) 
M°) J Pmc{k) 



+ 



R- 



Mo) 



d 3 kd 3 q(2ir) 



_ 3 Rs{k) (5(k) - Rs(J^ijk)l e ik.( q +(R z * Sz )( q )) k k 



P M c(k) 



J, VJ 



The integral over q in the last term above reduces to 5*(k). Therefore, using eq. (4.1) and 
the fact that by construction (8*6mc) = 0, we find that the ensemble averaged value of the 
last term is zero. Ensemble averaging the rest of the terms, we find that the Fisher matrix 
is given by 



1 



d 2 log C 



25 D (0) ds z ,i(k)ds* z>j (k) 



lJ \P L (k) + 



Rl 



FkqF x k 



e 



R 2 s(k) 
PAic(k) 



,(4.11) 



with x = q + (R z * s z )(q + q) — (R z * s z )(q) for an arbitrary q. For the second term above, 
we used both homogeneity and isotropy to simplify the ensemble average. The former brings 
out a Dirac delta function which cancels 1/5d(0) (and allows us to use an arbitrary q), while 
the latter tells us that the result must be proportional to 5ij. 

Combining the results above, we end up with the following iterative solution for s z : 



4(k) = s z (k)-W(k) 
R 2 z (k)P L (k) 



s(k) RiWPlW j. --i 



' . u Rs(k)(6(k) - Rs(k)di(k)) 

P M c(k) 



w-iw^i + ^^' ittJ 



m J xk 



2 R 2 s(k) 
" P M c(k) 



(4.12) 



where x = q+(R z *s z )(q + q)-(R z *s z )(q), x = q + (R z *s z )(q), and <5i = 5[q + (R z *s z )(q)] 
as before. The resulting expression above is written in a way suitable for implementation into 
a numerical code, since all integrals have been replaced by Fourier transforms. 5 To simplify 
our notation, we revert to denoting both Eulerian and Lagrangian wavevectors with k, unless 
stated otherwise. 

In Figure 1 we show the different quantities entering in eq. (4.12) for reference. In 
the actual numerical reconstruction we smooth out all wiggles in these quantities, so that 
wiggles are not introduced by hand, and only appear if present in the data. Note that the 
inverse Fisher matrix, W, is strongly weighted towards high k, which compensates the small-A: 
weighting coming from the Pmc m the denominator in (4.12). The Fisher matrix is obtained 
by performing the average in (4.12) over all N-body boxes described in Section 6, where for 
each box we have averaged over 100 randomly chosen q's. 

The transfer functions R z and R§ are 1 for small k, and deviate from 1 only for k > kj^L 
(see [7] for further discussion). The mode-coupling power, Pmc-, vanishes at small k, which 
is a direct result of the fact that 5\ and 6 are highly correlated for those scales [7] . 

4.4 Discussion 

The result of this section is equation (4.12), which gives a pseudo-optimal iterative solution 
for the Zel'dovich displacement field in Lagrangian space, given a non-linear density field in 
real space. There are, however, several important caveats to this result, which we discuss 
below. 

Note that (4.12) involves highly non-linear functionals of s z . Thus, one is not guaranteed 
that that iterative algorithm will find the global maximum of the Likelihood function, and 
not some local maximum. For example, the density field 5\ is insensitive to a swapping of 
pairs of particles. Such a swapping is in principle penalized by the term proportional to |s 2 | 2 
in C, (4.3), but it is not outright forbidden. Moreover, in the numerical calculation we do not 
want to introduce BAO wiggles by hand. Therefore, as mentioned before, we smooth out all 
wiggles from the quantities in (4.12), most notably - Pj_, and W. This automatically makes 
the algorithm not-completely-optimal (because of the modified W), and it does not converge 
to a true maximum of C (because of the modified Pl, entering in the first derivative of C). 



5 Note that the ensemble average in the Fisher matrix can be easily performed analytically, in the same way 
one calculates the analytical Zel'dovich power spectrum, Pz- However, unlike the Pz, the final Lagrangian 
space integrals in the Fisher matrix are not as trivial. So, we stick with an N-body approach to map k — > k 
through J'kqJ'^i^- One should note that whether or not one approaches the Fisher matrix by using N-body 
simulations, writing a code to do that mapping is still required for the first derivative of the Likelihood 
function. 
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Figure 1: We show the different quantities entering in the iterative solution for s z , eq. (4.12). 
The power spectra are divided by a smooth BBKS [9] power spectrum with shape parameter 
r = 0.15 in order to highlight the wiggles due to the BAO; and Pmc is multiplied by 5 for 
z = 1 to fit in the plot. Note that we have rescaled the inverse of the Fisher matrix W, to 
highlight the scale dependence. 



These same problems are further exacerbated because of our assumption of a Gaussian 5mc 
in writing down the Likelihood function. Even if the complications above were not present, 
one is always left with the problem of restoring the high-/c reconstructed power (see below), 
which has been killed by the Wiener-like filtering of the algorithm above. 

The above discussion implies that we should better start the algorithm, (4.12), as close 
to the true global maximum of C as possible. This is especially true, since we find the 
method to be very slowly converging. So, in Section 5 we construct a non-optimal (not even 
pseudo-optimal) but very efficient algorithm to provide us with a good initial guess for s z , 
which can then be fed into (4.12). With that guess for s z , we find that the algorithm above 
requires only 0(3) iterations to converge (see Appendix A for further details). 

Note that the discussion so far assumes that we are actually calculating the linear power 
spectrum directly from s z using (2.2). However, in the standard calculation of optimal 
quadratic estimators for Gaussian random fields, none of the problems above arise. Why 
is this the case? The reason is that for Gaussian fields one would integrate the likelihood 
function (4.3) over all reconstructed realizations, s z , and then optimize with respect to Pl, 
and not with respect to s z as we did above. That would result in a Pl, which is not affected 
by the smoothing applied by the Wiener-like filtering (e.g. (4.8)), because the Wiener filter 
drops out, as well as any possible bias introduced by a wrong initial guess for Pl, the R's, etc. 
However, integrating (4.3) with respect to s z is an immensely complicated task. Therefore, 
we have to make some approximation. And the approximation we do is finding Pl from 
the pseudo-optimal s z using (2.2). It is important, however, first to correct s z as much as 
possible at least for the smoothing introduced by the Wiener-like filtering (e.g. (4.8)). This 
is what we do next in an ad hoc manner. 

The discussion above implies that the residual between the obtained s' z and the true 
s z may still contain recoverable information about the BAO. In order to recover it, we write 



the relation between the true s z and the recovered s' z as: 

s z = R' z s' z + s' MC , (4.13) 

with {s' z s' MC ) = 0. A linear fit gives R' z = (s z ,is'* i )/(\s z \ 2 ). We find that indeed both R' z and 
s 'mc con tain important residual BAO information. So, we cannot use the trick of [3] to find 
the transfer function and the mode-coupling piece from other realizations (or slightly wrong 
cosmologies). Otherwise, we would be putting wiggles in our results by hand. Instead, we 
should recognize that if the cross-correlation between s z and s' z is very close to 1, this means 
that we have captured the complex phases of the perturbations correctly; and all that is left 
is to fix their amplitudes. Thus, we write our final result for the reconstructed s z as: 

s z (k) = A f (k)s' z (k) , (4.14) 

where Aj is an (inverse) Gaussian fit to A = ^/Pl/P' l , with 5d(0)P' l = \s' z \ 2 k. 2 . We use 
a Gaussian fit to A, since we recognize that A again contains wiggles, which we do not 
want to introduce by hand. When we perform the numerical calculations, we find that 
log(A f ) = (k/(1.19/i/Mpc)) 2 for z = and log(A f ) = (k/(0.97/i/Mpc)) 2 for z = 1, which are 
the best fit Gaussians for the scales relevant for the BAO: k < 0.4/i/Mpc. Note that ^4/(k) 
acts as the inverse of the Wiener-like filtering discussed above, thus restoring the high-k 
power in s z . In Section 6, we show results for the density field obtained by taking minus 
the divergence of s z , which is obtained by combining (4.12) and (4.14). But before we show 
those results, let us present the algorithm we use for initializing (4.12). 

5 Constructing a guess displacement field by modifying the standard re- 
construction scheme 

We find that the pseudo-optimal algorithm presented in the previous section converges very 
slowly - i.e. it is pseudo-optimal but is not efficient. Thus, we have to construct an efficient 
(although not even pseudo-optimal) algorithm to provide a guess displacement field for the 
pseudo-optimal algorithm of the previous section. Thus, the scheme presented in this section 
should be treated as providing a first-stage reconstruction, the results of which are then 
improved further by the second-stage reconstruction of the previous section. Inevitably, the 
scheme of this section is ad hoc, since we are only guided by the principle of choosing a 
scheme that is efficiently converging to a displacement field, which is well correlated with the 
true Zel'dovich displacement field. At the same time, we tried to make each of the steps of 
the algorithm below as reasonable as possible. 

Therefore, we will utilize the fact that there are several places where the standard 
reconstruction scheme can be improved. First, it would be good to write down the large- 
scale linear density field not through the Zel'dovich-like approximation given by eq. (2.4, 
2.6), but by using a more self-consistent approach to extract the large-scale linear density 
field. The other place where it can be improved is if we find a (more) optimal prescription 
for obtaining the filter Wq, which in the standard reconstruction is set to be a Gaussian with 
a smoothing scale given approximately by the rms particle displacements (~ 15 Mpc/h at 
z = 0). 

5.1 Obtaining a displacement field capturing the power in the non-linear density 

As we discussed in Section 3, the first stage reconstruction of this section is a backward 
scheme. In the spirit of a backward scheme, let us first find a displacement field d(q p ), which 
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nulls the power in the density, and therefore results in: 

S[x p - d(q p )} « , (5.1) 

where again x p and q p are the positions of the data particles in Eulerian and Lagrangian 
space, respectively. We will solve for the above equation iteratively below in a way analogous 
to [5] . The two main differences between our method and the iterative extension to standard 
reconstruction proposed in [5] is that: 1) we derive optimal filters (under a set of assumptions) 
with which we process the data instead of relying on filters chosen by hand; 2) we treat the 
large-scale modes consistently by not mixing Eulerian and Lagrangian coordinates. 

We know that most of the recoverable information for s z lies in 5± according to the 
model (4.1), and not in 5. Therefore as a first guess, we would like to match <5i: 

Si ~ S[q p + d{q p )\ as a first guess, (5-2) 

which requires obtaining 5\ from 5. We do that by noticing that the model (4.1) is linear. 
Thus, assuming that Smc is Gaussian, the optimal solution is to use a Wiener filter to obtain 
a recovered Si: 

S? c (k) = W s (k)S(k) , with W 5 (k) = ^SSfy , (5-3) 

where k is in Eulerian Fourier space. This Wiener filter can be cheaply obtained by calibration 
with simulations because of its small sample variance, in analogy with the transfer functions 
Rs and R z discussed above. We plot the filter Ws in Fig. 2 for redshifts z = and 1. Note 
that the filter decays at about the non-linear scale (given by kpjL = 0.25/i/Mpc for z = and 
kNL = 0.74/i/Mpc for z = 1) as one can expect. It also features slight wiggles corresponding 
to the BAO wiggles. In order not to introduce those wiggles by hand, we smooth them out 
by fitting a B-spline function to Ws, and use that result in our numerical calculations. 

To obtain d(q p ), we apply the iterative method described in [5], smoothing intermediate 
results for the density with Ws (see Appendix B for further details). In this iterative scheme 
we start with the true particle positions, which correspond to S. We smooth S with Ws to 
obtain a reconstructed Si, and then derive the displacement field corresponding to it (using 
(2.2)). That displacement field approximately matches d in (5.2). We then displace the 
particles by minus that displacement field; recompute the density; smooth it with Ws; derive 
a new displacement field; and move the particles again. We iterate this procedure N times, 
where fiducially N = 5. After these 5 iterations, the power spectrum of the density field 
corresponding to the newly obtained particle positions (x p ) is only about 10% (for z = 0) 
of the power in Sl at k ~ 0.75/i/Mpc and is strongly suppressed at lower k. Thus, the 
newly obtained particle positions, x p , result in a more or less uniform density field (up to 
k ~ 0.75/i/Mpc). So, x p should be close to the Lagrangian position of the particle p apart 
from some short-scale corrections, which we will quantify later on (see eq. (5.4)). 

That implies that for the purposes of the acoustic peak reconstruction, we satisfy 
eq. (5.1), where d(q p ) is given by d p = x p — x p , which is the total displacement of 

particle p after N iterations. Therefore, the pairs (x p ,d p ) equal (q p ,d(q p )) up to some 
short-scale corrections. So, all we have to do next to obtain d(q) is to interpolate those pairs 
to a grid. This can be done in variety of ways. Our method for accomplishing that is given 
in Appendix C. It is an iterative method (the fiducial number of iterations (M) we choose 
for the method is M = 5), which produces d(q) on a grid in Lagrangian space. 
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z = 1 




k[h/Mpc] 
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Figure 2: The density and displacement Wiener niters discussed in the text. The filter, 
Wd, is calculated for N = M = 5. In order to remove the BAO wiggles from Wjj, in 
our reconstruction scheme we use a Gaussian fit to Wd, denoted by Wb,exp- Note that k 
corresponds to Eulerian Fourier space for Wg, and to Lagrangian Fourier space for Wd and 

W D ,exp- 



5.2 Extracting the Zel'dovich displacement field 

Note that in the steps above we obtained some d satisfying (5.1), starting from a guess satis- 
fying (5.2). So, next we need to relate s z to d, compensating for any short-scale corrections, 
and non-linearities. We can always write that relation as follows (in Lagrangian Fourier 
space) : 



d(k) = R D (k)s z (k) + smc(K) 



(5.4) 



with {sz^smcj) = 0. That allows for a linear fit for the new transfer function, Rd, giv- 
ing Rd = (s z • d)/(|s z | 2 ), which can be calibrated after applying the algorithm above to 
simulations. The equation above in principle should capture any problems arising from the 
inversion of (5.1), as well as any non-linear effects and noise. To make progress however, we 
assume that the "mode-coupling" term, smc, is Gaussian. This tells us that the optimal s z 
can be obtained through Wiener filtering: 



a£(k) = W D (k)d(k) , with W D (k) 



(^(k)-d(k)) 
(|d(k)P) 



(5.5) 



We put a superscript A to highlight that the above equation recovers only the low-fc modes. 
Note that the optimal Wd depends on the number of iterations A^ and M, and as in the 
case of W$ can be cheaply extracted from simulations. Using equation (2.2), the low-fc linear 
density field is given by 



'(k) 



-is, 



-iWo(k)d(k) 



(5.6) 



where "new" denotes that the linear density reconstructed with the scheme proposed in 
this paper. Note that the resulting large-scale reconstructed linear density field, 5^ w , is 
consistently derived in Lagrangian space. This is in contrast to standard reconstruction, 
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where the large-scale modes are reconstructed using a Zel'dovich-like approximation (2.4, 
2.6), which mixes Lagrangian and Eulerian coordinates. 

We plot the filter Wd (see next Section for description of the N-body simulations we 
use) in Fig. 2 for z = and 1, resulting after applying the scheme above with N = M = 5. 
Similar to W$, Wd decays roughly on the non- linear scale, although the decay is steeper than 
for Ws- The filter features strong wiggles corresponding to the BAO wiggles. In order not to 
introduce those wiggles by hand into the reconstruction scheme, we remove them by fitting 
a Gaussian to Wd, and using that fit instead of Wd in our calculations (e.g. in eq. (5.6)). 
We show the Gaussian fit to Wd in Fig. 2 for reference (the curves labeled WD,exp)- The 
Gaussian scale in Woexp is 2.80Mpc//i for z = and 1.54Mpc//i for z = 1. 

5.2.1 Recovering the short-scale modes 

The procedure which we used above to obtain s^ is definitely non-optimal (To remind the 
reader the pseudo-optimal scheme is given in eq. (4.12)). Thus, the residual d — s^ contains 
important information about s z . There are many choices of how one can recover the short- 
scale modes. We find that a method that works well is to move a uniform field of particles 
(with positions q u ) with the short-scale residual displacements, which in Fourier space are 
given by d s (k) = (1 — W,o(k))d(k). Thus, we have: 

C W = %* + d s (q u )} , (5.7) 

which corresponds to the density field obtained by pulling back the data particles by the 
large-scale s^ as done in standard reconstruction. Note that all of the relevant information 
has been transformed (although arguably this is a lossy transformation) from 5 to d, and 
therefore we do not need to add the residual 5 R from (B.2). Indeed adding it, only degrades 

the performance of our method, while increasing N makes S R from App. B increasingly 
irrelevant. 

5.2.2 Putting it together 

The final linear density approximator for this first stage of the reconstruction is then given 
by 

5l™(q) = 5l™(q) + 5*™( q ). (5.8) 

Combining the above equation with (2.2), we can find the displacement field s° cw correspond- 
ing to 5 L cw (q): 

k 

„new ■ rnew /r n"\ 

S z =l ^2°L ■ (5-9) 

It is this s^ ew , which we then use as an initial guess for the pseudo-optimal second stage of 
the reconstruction, given by (4.12). 

6 Numerical results 

To gauge the performance of our reconstruction scheme against standard reconstruction, we 
ran 26 N-body simulations using the GADGET-2 code [10], with 256 3 particles each, and a 
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Figure 3: Cross-correlation function between the linear density field and the density field 
labeled in the legend. The "NL" line corresponds to the cross-correlation coefficient between 
the linear and non-linear densities; while the "Rec" curve corresponds to the reconstruction 
scheme proposed in this paper. The curves labeled with different values of E, correspond to 
the standard reconstruction result with a Gaussian smoothing on the scale E. 

box size of 500Mpc//i. The cosmological parameters we use are (in the standard notation) 
as follows 



(Pb, ^matter, «a, h, n s , <t%) = (0.046, 0.28, 0.72, 0.70, 0.96, 0.82) 



(6.1) 



The initial conditions are set at a redshift of 49, using the 2LPT code provided by [11]. 

Below we show results from the two-stage reconstruction scheme presented in this paper. 
The results for the first stage of the reconstruction are denoted by "IRec" in the figures that 
follow. That corresponds to the reconstructed linear density field from Section 5 (given by 
eq. (5.8)). That reconstructed linear density field is used to construct a guess displacement 
field s" ew , (5.9), which is then fed into the iterative solution for s z given by (4.12) and then 
into (4.14). The final reconstructed Zel'dovich displacement field, s z , obtained from (4.14) 
is converted to a reconstructed linear density field, again using (2.2). That result is denoted 
by "Rec" in the figures that follow. We refer the reader to Appendix A for further details on 
our numerical implementation of (4.12). 

In Fig. 3 we show the cross-correlation between the linear, Sl, and the non- linear density 
5 (curve denoted by "NL"), as well as the one between the linear and reconstructed densities. 
The curves for the standard reconstruction result are denoted by the value of the smoothing 
length E entering in Wq- The cross-correlation between our reconstructed linear density fields 
and the true linear density field is shown with the curves denoted by "IRec" and "Rec". 

One can see that standard reconstruction provides markedly different results for different 
E. In the limit of very large E, Wg ~^ and we obtain 5\ = and S s = 5. This implies that 
the cross correlation between 5^ and 5l approaches the one between the linear and non-linear 
density fields - a behavior which one can clearly see in Fig. 3. In the limit of very small E, 
we obtain 5 S = and a 5\ which according to the discussion at the end of Section 2 can be 
roughly approximated by the density in the ZA. Therefore, the cross correlation between 6^ 
and Si can be approximated by the one between the linear and ZA density, which is in turn 
well approximated by the cross-correlation between the linear and non-linear density fields 
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Figure 4: Shown are the linear and non-linear (obtained from N-body simulations) density 
power spectra, as well as the power spectra of the reconstructed density fields. The standard 
reconstruction result is given by the curve denoted by E; while the result of this paper 
is denoted by "Rec". The errorbars show the 2-sigma errors in the average of the linear 
power-spectrum; and are representative of the errors in the rest of the power-spectra. The 
/c-positions of the errorbars indicate the binning of the power spectra, while the curves were 
obtained after splining the result. Note that ideally the reconstructed density power spectrum 
should match the linear power, i.e. should not be simply within the errorbars. In practice, 
however, the BAO peak is determined by the period of the oscillations and is not affected by 
the broader-band features. 



(e.g. [3]). This behavior can be seen in the appearance of the step at low k observed for 
E = 3Mpc//i in Fig. 3. Therefore, by varying E we can trade between reconstructing either 
the large-scale modes (k < 0.2/i/Mpc for z = 0); or the small-scale modes (k > 0.2/i/Mpc for 
z = 0) at the expense of the large-scale modes. 

Our reconstruction scheme does not treat the large-scale modes in the Zeldovich-like 
approximation of eq. (2.6), and instead tries to reconstruct the true linear density through 
eq. (2.2). Thus, from Fig. 3 one can see that both stages of our reconstruction scheme recover 
both the large-scale modes along with the short-scale modes with high fidelity as measured 
by the cross-correlation. The first stage of our reconstruction is analogous to the standard 
reconstruction scheme, apart from our treatment of the large-scale modes. This implies that 
the main gain of our proposed methods comes from the fact that we obtain the large-scale 
linear power by applying eq. (2.2), i.e. by explicitly avoiding the mixing up of Eulerian and 
Lagrangian coordinates done in the standard scheme. 

Note that the high-E curves have slightly better cross-correlation at low k than our 
reconstruction scheme. However, this does not help the reconstruction of the acoustic peak 
(as we will see below), because it is at the expense of a much poorer reconstruction of the 
small-scale modes. 
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Figure 5: Shown are the 2-pt matter correlation functions corresponding to the power 
spectra in Fig. 4. Due to the missing large-scale power due to our relatively small box sizes, 
the 2-pt functions are shifted down relative to the true 2-pt functions, and a slight extra 
slope is introduced. However, since we compare results obtained from the same N-body 
simulations, these large-scale effects should not affect our conclusions. To reduce sample 
variance, the z = plot is obtained after binning the data in radial bins of width 3.9Mpc//i; 
while the bins for z = 1 are of width 1.95Mpc//i (apart from the NL curve which is binned 
in twice larger bins). 

In Fig. 4 we show the linear and non-linear density power spectra, as well as the power 
spectra of the reconstructed density fields. The corresponding 2-point functions are shown 
in Fig. 5. From Fig. 3 we chose £ = 10Mpc//i for z = 0, and £ = 5Mpc//i for z = 1 as 
fiducial values for the smoothing length in the standard reconstruction scheme. Note that 
the reconstructed 2-pt function is rather insensitive to the choice of £ - the fiducial £ values 
corresponding approximately to the best-reconstructed acoustic peak with that scheme. 

Ideally, the reconstructed density power spectrum should match the linear power spec- 
trum. The standard reconstruction scheme recovers the BAO wiggles up to 6 k ~ k^L = 
0.25/i/Mpc at z = 0. The reconstruction scheme proposed in this paper recovers the baryon 
wiggles up to k rs 0.4/i/Mpc at z = after the pseudo-optimal second stage of the recon- 
struction (and up to k ~ 0.35/i/Mpc after the first stage), which corresponds to roughly two 
more well-reconstructed wiggles in the power spectrum. 7 Therefore, the peak in the corre- 
sponding 2-pt function (see Fig. 5) shows marked improvement compared with the result 
from standard reconstruction. 

The improvements for z = 1 in the 2-pt function and the power spectrum with our 
scheme are not that stark when compared to the results from the standard reconstruction 



6 This maximum k for which we see wiggles in the power spectrum depends on the smoothing scale E. One 
can reduce E by about a factor of 2 and recover wiggles with standard reconstruction up to k ~ 0.35/i/Mpc 
for z = 0. However, that is at the expense of reducing the power for scales k ~ 0.1/i/Mpc, due to the fact 
that standard reconstruction treats the large-scale modes in a Zel'dovich-like approximation. This results in 
an acoustic peak, which is more smeared than for our fiducial choice of E. 

7 One might worry that the first stage of our reconstruction gives better match to the linear power spectrum 
at k ~ 0.15/i/Mpc than our second stage. However, this seems accidental, since if one "fixes" the large-scale 
power in the first-stage result (e.g. by multiplying by a broad-band Gaussian decay), this result would go 
away. 
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scheme. Still, one should note that the power spectrum that we recover for z = 1 lies very 
close (within a couple of percent for our choice of binning) to the true linear power up to 
k m 0.4/i/Mpc (after the second stage of our reconstruction scheme), testifying to the quality 
of the reconstruction scheme proposed in this paper. 

7 Summary 

In this paper we derive a pseudo-optimal scheme for reconstructing the baryon oscillations 
in the matter power spectrum. The proposed scheme has no free parameters (apart from the 
grid size used for the Cloud-in-Cell assignments). It is based on a simple model (see (4.1)) 
for the non-linear density proposed in [7]. 

The reconstruction scheme of this paper is made up of two stages. The first consists 
in finding a good guess for the Zel'dovich displacement field. That field is then fed to the 
second stage, which solves for the Zel'dovich displacement field that maximizes the Likeli- 
hood function for our model for the non-linear density field. The first stage requires O(10) 
iterations, while the second stage - 0(3). 

The main differences between the reconstruction scheme proposed here and the standard 
reconstruction method is 1) in using optimal filters in both stages of the reconstruction; 
and 2) in our treatment of the large-scale modes. The first point implies that our method 
has no free parameters such as the smoothing scale entering in the standard reconstruction 
scheme. Regarding the second point, standard reconstruction obtains a large-scale density 
field which is in fact evaluated using a Zel'dovich-like approximation. That choice degrades 
the reconstruction of the large scales especially for the low values of the filtering scale which 
are needed if one wants to recover the short scales better. Such a trade-off is in practice 
absent in our method, which also avoids the mixing up of Lagrangian and Eulerian space 
coordinates, which is done in the standard scheme. 

We find that for z = the reconstruction scheme proposed in this paper recovers 
the baryon wiggles up to k ~ 0.4/i/Mpc, compared with k ~ 0.25/i/Mpc for the standard 
reconstruction scheme. This implies that we recover roughly two more wiggles in the power 
spectrum, which results in a markedly improved sharpening in the reconstructed acoustic 
peak compared with the result from standard reconstruction. However, we postpone a proper 
error analysis for the reconstructed acoustic peak position to future work. 

We find that even the first stage of our scheme yields excellent reconstruction, recovering 
the baryon wiggles up to k ~ 0.35/i/Mpc for z = 0. Whether this result will persist for biased 
tracers, remains to be seen. 

One can envision numerous possible improvements of the scheme proposed in this paper. 
One such improvement is using the model in [7] for the non-linear density field based on 
second-order Lagrangian perturbation theory, instead of the one used in this work based on 
the ZA. We postpone such analysis, as well as the application of the proposed method to 
biased tracers in redshift space, to future work. 

A Implementing (4.12) numerically 

To use (4.12), we need to map x to q between the Fourier transforms F^qF^Qfo), where 
Q can be read off from (4.12). This is done by displacing a uniform grid of particles at q u 
by (R z * s z )(q u ) (or by (R z * s z )(q u + q) — (R z * s z )(q) for the Fisher matrix). Then we 
interpolate the quantity F~ k Q(k) to the particle positions using three-linear interpolation, 
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and assign the obtained values to the initial q u . The resulting field corresponds to evaluating 
J r a T fc Q(fc) in Lagrangian space and is therefore ready for feeding into J 7 ^. 

Another thing to note is that to make the algorithm stable we truncated Q in the first 
derivative of C (the quantity in the inner square brackets in the first line of (4.12)) at 1/2 the 
Nyquist wavevector, corresponding to k ~ 0.8/i/Mpc. This is reasonable since that quantity 
involves derivatives in Eulerian space. We checked the sensitivity of our method to this 
choice of cutoff, by pushing the cutoff to 1/i/Mpc. This resulted in only a minor change of 
the amplitude of the reconstructed Pl at the high k end, which did not affect the 2-pt or the 
cross-correlation functions. 

Note that for all iterations of (4.12) we evaluate W at the maximum of the likelihood 
function, using the result shown in Fig. 1, but with the wiggles removed. Given the excellent 
guess coming from the first stage of our reconstruction, this does not affect the final result, 
at the same time speeding up the calculation appreciably. 

One final modification that we needed to do to (4.12) to make it more stable, was to 
multiply W by a constant, call it 7. This is a standard trick used with the Newton- Raphson 
algorithm when applied to non-quadratic functions 8 , such as our log£. We find that an 
initial 7 = 0.2 works well for our guess s z coming from the first stage of the reconstruction, 
described in Section 5. We restart an iteration with a 7 reduced by a half whenever the ratio 
of the rms displacement corresponding to s z to the rms displacement corresponding to the 
correction (the term proportional to W in (4.12) but before multiplying it by 7) is smaller 
than 0.95 times the ratio of the previous iteration. 

We find that this ratio of rms displacements is always bigger than one and quickly 
improves. We stop the iterations in (4.12) when that ratio exceeds 10 for z = and 30 for 
2 = 1, since we find that further iterations do not improve the result further. This stop 
condition is reached after 0(3) iterations, given our initial guess for s z . 

B Solving eq. (5.1) for the displacement field starting with eq. (5.2). 

In this section, we solve eq. (5.1) for d{q p ) starting with eq. (5.2). First, we construct a 
displacement field on a grid, which is initially set to zero: g(°' = 0. We then apply the 
following iterative scheme on g in a way analogous to the method described in [5]: 



g^(k)=i-^W s (k)5^- 1 \k), (B.l) 



where 

4 n) = 6[xp>] , 4 n) = x p n ^ - g {n \xf-^) . (B.2) 

In writing down the equations above we applied the Wiener filter, Wg, entering in eq. (5.3). 
We evaluate g^ n ' at x p n by using three- linear interpolation. The initialization is given by 
x p = x p , i.e. we start with the true particle positions, which implies that 5 R =5, and that 
g' '(k) approximately matches d in (5.2). 

As discussed before (after eq. (5.3)), each successive iteration in the above algorithm 
kills the large-scale power in 5r to higher and higher k. So, after 5 iterations, the power in 5r 
for z = is only about 10% of the power in 8l at k ~ 0.75/i/Mpc and is strongly suppressed 



In fact, the presence of 7 may be required even for quadratic functions to make the algorithm convergent, 
since we use the expectation value of the Hessian, and not the Hessian itself. 
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at lower k, which implies that for the purposes of the acoustic peak reconstruction, we satisfy 
eq. (5.1). 

The total displacement by which we have moved particle p after iV iterations is given 
by: 



N 



4*0 = £^(^-1)), (B . 3) 



n=l 



which is a vector which should be anchored to the particle position after the iV-th iteration, 

Usp 

C Interpolating the displacements to a grid 

We would like to construct an approximation for d(q) starting from the pairs (x p ,d p ') 
which approximately equal (q p ,d(q p )) up to some short-scale corrections. One can do that 
in a variety of ways, e.g. using Voronoi tessellation to perform the interpolation to a grid. 
For the illustrative purposes of this paper, however, we choose a simple iterative algorithm: 
We use a weighted CiC assignment to construct the mass-weighted displacement field on a 
grid. Dividing that by (1 + 5 R ) we obtain a guess displacement field on a grid. Call that 
G^'(q). Then we perform the following iterative procedure: 

dp N - m ) = d p N '-- 1 ) - CW^) , W ith df.°) = 4*0 , (C.l) 

where G^ m ' is obtained by dividing the mass- weighted displacement field corresponding to the 
pairs (x p ,d p ' m ) by the density (l+<^ ). To avoid singularities, we set the displacement 
in zero density regions to zero. After M iterations (in our case we find that M = 5 works 
well enough), the approximate displacement field in Lagrangian space is given by 

M 

d{q)=Y,G {m) {Q)- (C2) 

m=l 
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